Monte Carlo study of the two— dimensional site— diluted dipolar Ising model 
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By tempered Monte Carlo simulations, we study 2D site-diluted dipolar Ising systems. Dipoles 
are randomly placed on a fraction x of all sites in a square lattice, and point along a common 
crystalline axis. For Xc < x < 1, where Xc ~ 0.79(5), we find an antiferromagnetic phase below a 
temperature which vanishes as x ^ Xc from above. At lower values of x, we study (i) distributions 
of the spin-glass (SG) overlap q, (ii) their relative mean square deviation and kurtosis and (Hi) 

/-£/, where is a SG correlation length. From their variation with temperature and system size, 
we find that the paramagnetic phase covers the entire T > range. Our results enable us to obtain 
an estimate of the critical exponent associated to the correlation length at T = 0, l/v = 0.35(10). 
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I. INTRODUCTION 

In the last years, there has been a renewed interest in 
systems of interacting dipoles (SIDs). This is in part due 
to recent advances in nanoscience^ which make realiza- 
tions of assemblies of magnetic nanoparticles available)^ 
Empirically, these systems show a rich collective behav- 
ior in which the dipole-dipole interaction plays a key role 
that can be observed at low (but experimentally acces- 
sible) temperatures. Dipoles forming crystalline arrays 
exhibit long-range ferro or antiferromagnetic order that 
depends crucially on lattice geometry^ because of geo- 
metric frustration caused by the spatial variations of the 
directions of dipolar fields. Two-dimensional (2D) arrays 
of cobalt-ferrite and Co nanoparticles placed on hexag- 
onal arrays have been found to exhibit in-plane short- 
range ferromagnetic order.— On the contrary, arrays on 
a square lattice composed of MnAs ferromagnetic nan- 
odisks epitaxially grown on a substrate exhibit collinear 
AF patterns.— 

Magnetic ordering of SIDs depends also on anisotropy. 
On the one hand, dipolar-dipolar interactions create ef- 
fective anisotropics that in square lattices, for example, 
push spins to lie on the plane of the latticed On the 
other hand, magnetocrystalline site-anisotropy energies 
of the crystallites that form the nanoparticles are often 
greater than dipolar-dipolar interparticlc energies. This 
is the case of the arrays of MnAs ferromagnetic nan- 
odisks we mention above, that behave as a system of Ising 
dipoles with their magnetic moment rigidly aligned along 
the in-plane crystalline easy axes of the nanodisks.^ In 
such a case, the resulting magnetic order depends on the 
competition of dipolar and anisotropic energies. Crys- 
talline Ising dipolar systems (IDSs) are reasonable mod- 
els for these planar systemsi^ Some ferroelectricsr^ insu- 
lating magnetic salts as LiHoF4, as well as some three- 
dimensional (3D) crystals of organometallic molecules^i 
are known to be well described by arrays of IDSs.— 

SIDs in disordered spatial arrangements are particu- 
larly interesting. The presence of spatial disorder, to- 



gether with the geometric frustration generated by dipo- 
lar interactions, gives rise to random frustration that may 
result in SG behavior. In fact, some non-equilibrium SG 
behavior (like time dependent susceptibilities and mem- 
ory effects) has been observed in experiments with sys- 
tems of randomly placed nanoparticles or very diluted 
magnetic crystalsJ^ii^ Furthermore, Monte Carlo (MC) 
simulations have given clear evidence of the existence of 
a transition at finite temperature Tsa from a paramag- 
netic to an equilibrium SG phase in systems of randomly 
oriented axis dipoles (RADs) placed either on fully oc- 
cupied or on diluted simple cubic (SC) 3D lattices, and 
TsG = instead for 2D square latticesJ^ Recent numer- 
ical work has reported a SG transition in a model of par- 
allel axis dipoles (PADs) placed on a lattice that approx- 
imates that of the dilutedi^ Li IIo2;Yi_^ F4, a material 
for which such a transition has been reportcdfi^ (albeit 
not without some controversyi^) . By MC simulation the 
whole phase diagram of site-diluted PADs placed on a 
3D SC lattice has been obtainedi^ as a function of the 
concentration x. It includes a SG phase for < x < 0.65 
which, strikingly, has been found to behave marginally, 
that is, it has quasi-long range order, as in the 2D XY 
model This is contrary to theoretical expectationsr^i^i 
that SG systems with long-range interactions may be- 
have as short-range Edwards-Anderson (EA) models;^ 
which in 3D are believed to have a SG phase with a non- 
vanishing order parameter (according to the RSB^ or 
droplet^i. pictures of SGs). 

Then, in order to get a deeper understanding of SG 
systems beyond the already extensively studied random- 
bond models with short-range interactions, it makes 
sense to analyze the behavior of the 2D PAD model and 
compare it with the short-range 2D EA model. The lat- 
ter had been found to have an algebraic divergence at 
TsG = with critical exponent^^ l/v = 0.50(5), al- 
though more recent simulations for larger systems and 
lower temperatures give a value of 1/v = 0.29(4) for 
Gaussian interactionsi^ 

Our purpose is the study by MC simulations the phase 
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diagram of a site-diluted system of magnetic dipolcs. 
They are placed at random on the sites of a square lat- 
tice and point up or down along a given principal axis. 
Since in the limit of low concentrations every detail of 
the lattice is expected to become irrelevantji^ our results 
have direct connection with some of the work we describe 
above. Our intention is to search for the temperature Tsg 
of a possible SG transition and study the related diver- 
gence of the correlation length. Further, we aim to study 
whether the diluted PAD model belongs to the same uni- 
versality class recently conjectured?^ though not reliably 
shown by MC simulations?^ for the set of 2D EA Ising 
models with varying quenched disorder. 

The plan of the paper is as follows. In Section|II]we de- 
fine the model and give details on the parallel tempered 
Monte Carlo (TMC) algorithm^S used for updating. We 
also define the quantities we calculate. They include the 
spin overlap^ q and a correlation lengt h^^i^^ . In Sec- 
tion nil Al we give results for the dipolar AF phase for 
X > Xc, where Xc = 0.79(5), as well as for its nature 
and boundary. In Section IIII Bl numerical results are 
shown for distributions of q and / L at x — 0.2 and 0.5. 
We examine the evidence against the existence of a finite 
temperature SG phase transition when x < x^'- (i) the 
mean values (| q |) and (g^) decrease faster than alge- 
braically with L as L increases for T/x > 0.3, (ii) double 
peaked, but wide, distributions of q/ {\ q \) change with L 
for temperatures as low as T/x — 0.4, and (in) kurtosis 
and £,l/L decrease with L at all T and do not cross, as 
it would be expected for a finite temperature transition. 
Scaling plots for g and S,l/L arc given in Section BlI CI 
Our results are consistent with a ratio ^l/L that diverges 
with exponent 1/u = 0.35(10). Results are summarized 
in Section ITVl 



II. MODEL, METHOD, AND MEASURED 
QUANTITIES 

A. Model 

We treat site-diluted systems of Ising magnetic dipoles 
(also named spins in this paper) on a 2D square lattice. 
At each lattice site a dipole is placed with probability x. 
Then, the number of spins on the lattice is less than 
[L is the lateral size of the lattice) approximately 
by a factor x. Site i is said occupied if it contains one 
spin. All dipoles are parallel and point along the Y axis 
of the lattice. This axis shall be called spin axis. The 
Hamiltonian is given by, 

H = -^T^jai(Tj , (1) 

where the sum runs over all occupied sites i and j except 
* = = il on ^iiy occupied site i and 

T., =£.(a/r,,)'(l-34/r2.). (2) 



If Yij is the vector joining sites i and j, then = \\Tij\\ is 
its modulus and yij its Y component. Ea is an energy and 
a is the lattice spacing. In the following all temperatures 
and energies shall be given respectively in units of Ea/ks 
(ks is the Boltzmann constant) and Sa- 

Due to the long-range nature of the dipolar interac- 
tions, we are able to simulate on rather small lattice sizes 
(L < 32). 

Strength Tij is the usual long-range dipole-dipole in- 
teraction. Note that Tij signs are not distributed at ran- 
dom, but depend on the orientation of vectors on the 
lattice. Randomness in our model arises only through the 
introduction of the probability x for placing dipoles. This 
is to be contrasted with random-bond EA Ising models 
with bond strengths Jij = ^ij/^ij where Sij are chosen 
at random from a bimodal or Gaussian distribution with 
zero mean^i and /i is a real exponent. This is why PADs 
exhibit AF order at high concentration in contrast to 
these models that do not. Similar statements apply when 
our PAD model is compared with a random-axes dipolar 
model (RAD), in which Ising dipoles lie along directions 
chosen at random for each sitCf^ 



B. Method 

Periodic boundary conditions (PBC) are imposed. 
Spins on occupied sites i have been allowed to interact 
only with spins j within an L x L squared box centered 
on site i. This method unambiguously defines the vector 
Vij to be used in ^ and also excludes interactions with 
spins belonging to the repeated copies of the lattice that 
appear beyond the boundary. Because of the long-range 
nature of dipolar interactions, contributions from beyond 
this box would have been taken into account (for exam- 
ple by means of Ewald's summations^) if spins were to 
form ferromagnetic domains. They do not do so in our 
PAD model. In all simulations presented in this work we 
have found Txf ^ 1, where xf is the ferromagnetic sus- 
ceptibility. Therefore those contributions do not affect 
the thermodynamic limit regardless of the kind of phase 
the system lies (paramagnetic, AF or SG). Some details 
on this point are found inJ^ 

In order to circumvent large energy barriers that could 
slow down the evolution of the system, in particular 
from certain states representing minima of the energy 
(mainly at low temperatures), we have used the TMC 
algorithm.— It consists in running in parallel a set of 
n identical systems at equally spaced temperatures Ti, 
given by T,=Ti-{i- 1)AT (i = 1, • • • , n and AT > 0) 
where each system i is cyclically allowed to exchange its 
state with system i + 1. Each system evolves indepen- 
dently by use of the standard single-spin-flip Metropolis 
algorithm^ and whenever a single flip is accepted, all 
dipolar fields throughout the entire lattice are updated. 

In detail the procedure is as follows J^ii^ (V ^ cycle on 
i is run from i = 1 to i = n; (2) when the cycle arrives at 
system i, 8 Metropolis steps are applied on it; (3) next. 
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TABLE I: Simulation parameters, x is the probability for 
sites to be occupied with a magnetic dipole; L is the lateral 
lattice size; AT is the temperature step in the TMC runs; Ti 
and Tn are the highest and lowest temperatures, respectively; 
Nr is the number of pairs of quenched disordered samples; to 
is the number of MC sweeps. The measuring time interval is 
\to, 2to] in all cases. 



X = 0.2, AT = 0.02, Ti = 0.6, to = 4 x 10^ 


L 


8 12 16 20 24 32 


Tn 


0.04 0.04 0.04 0.04 0.04 0.08 


Nr 


2400 550 1500 650 700 200 


X = 0.5, AT = 0.05, Ti = 2, T„ = 0.1, to = 8 x 10'' 


L 


8 16 20 24 


Nr 


2500 2500 350 250 


X = 0.6, Tn = 0.2, to = 4 X 10** 


L 


8 16 20 24 


Ti 


3 3 2 2 


AT 


0.2 0.2 0.1 0.1 


Nr 


1200 300 300 300 


X = 0.7, AT = 0.1, Ti = 2, T„ = 0.2 


L 


8 16 20 24 


to 


8 X 10^ 8 X 10^ 4 X 10"3 4 x 10*' 


Nr 


4200 2200 400 100 


X = 0.8, AT = 0.1, Ti = 3, T„ = 0.2 


L 


8 16 20 24 


to 


8 X 10^ 8 X 10^ 4 X 10^ 4 x 10'' 


Nr 


4500 1200 500 350 


X = 0.86, AT = 0.1, Ti = 3, T„ = 0.2, to = 8 x 10^ 


L 


8 16 20 24 


Nr 


3000 400 300 70 


X = 0.9, AT = 0.1, Ti = 3, T„ = 0.2, to = 8 x 10'^ 


L 


8 16 20 24 


Nr 


2000 250 250 800 



a chance is given to systems i and i + 1 to exchange their 
configurations (note that at this moment system i + 1 
has undergone 8 Metropohs steps less than system i). 
The exchange is accepted with probability Ptmc = 1 if 
SE = Ei-Ei+i < or Ptmc = cxp{-Af]6E) otherwise. 
Here is the numerical value of Hamiltonian ([1]) for 
system i and A/3 = l/T^+i — l/T^; (4) 8 Metropohs steps 
are apphed on system i + 1 (regardless of the fact that 
the previous exchange have or have not been performed); 
(5) the above exchange is tried between systems i + 1 and 
1 + 2; (6) the cycle ends after the 8 Monte Carlo steps for 
i = n, after which no configuration exchange is tried. 

Sincei^ in 3D Tsg ^ x and the purpose of TMC is to 
overcome energy barriers that could be as high as Tsg, 
then we found necessary to choose the highest temper- 
ature Ti > 2x. It is also important to take AT small 
enough to allow frequent state exchanges between sys- 
tems. This is fulfilled by taking AT < T/VcIlV where 
Cs is the specific heat per spin. We choose appropriate 
values for AT from inspection of plots (not shown) of 
the specific heat vs T in preliminary simulations of the 
smaller systems 



Initially the n configurations were completely disor- 
dered. For details on how we chose equilibration times Iq 
see Section III CI Time to is particularly large outside the 
AF region, varying from at least 4 x 10^ MC sweeps for 
X = 0.7 and a number of dipoles N > 300 up to 4 x 10'' 
sweeps for x ~ 0.2 and N — 200. Instead, to in the AF 
zone is as low as 8 x 10^ for x > 0.86. Thermal averages 
were calculated over the time range [to,2to]- We further 
averaged over A^^ samples with different reahzations of 
disorder. Each realization was run twice to permit the 
calculation of overlapping parameters (see Section III Cp . 
Values of the parameters for all TMC runs are given in 
Table I. 



C. Measured quantities 

Measurements were performed after two averagings: 
first over thermalized configurations and secondly over 
different realizations of the quenched disorder. 

We begin by presenting the specific heat. It was ex- 
tracted from the slope of the energy as a function of the 
temperature. 

As for the staggered magnetization, also for a PAD 
model on the square lattice we find appropriate to define 
it as^ 

m = N-'Y.i-ir''^' (3) 

i 

where Xi is the X coordinate of site i. We calculated 
the probability distribution P{m), as well as the mo- 
ments mi = (|m|), 7712 = (m^), and mi = {m'^), where 
(...) stands for the above-defined double averages. From 
these moments we calculated the kurtosis (known also as 
Binder's cumulant) of P{m) as = {'i — m4^/m\)/2. All 
these quantities have proven to be good signatures for 
possible AF-paramagnetic phase transitions. 

In order to look for SG behavior, we also calculated the 
Edwards- Anderson overlap parameter between two inde- 
pendent equilibrium configurations obtained from a pair 
of identical replicas evolving independently in time^^ 

g = A^-i^</>,, (4) 

3 

where 

0. = ) , (5) 

a^p and (7^^ being the spins on site j of replicas labelled 
as (1) and (2). Clearly, g is a measure of the spin config- 
uration overlap between the two replicas. As done for m, 
we also calculated the probability distribution P{q) as 
well as the moments qi = {\q\), q2 = {q'^), and q4 = (g^). 
The SG susceptibility XSG is given by Nq2. Finally, we 
also make use of the relative mean square deviation of q, 
= 92/?! — 1, and kurtosis g = {3 — g4/q'|)/2. 
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FIG. 1: (Color online) Semilog plots of 52(^0, i) and q2 vs time 
t (in MC sweeps) for systems of 24 x 24 sites and concentration 
X = 0.5 at the values of T shown in the legend. Here, 52 comes 
from averages over time, starting from an initial random spin 
configuration at t = 0. Here to = S x 10® MC sweeps. Data 
points at time t from an average over the time interval [t, 1.2t] 
and over 500 system samples. 



FIG. 2: (Color online) Phase diagram of the 2D PAD model, 
o stand for the Neel temperature Tjv and x for temperatures 
below which we cannot completely rule out a SG phase (see 
Section HIH) . The full line for the phase boundary between the 
paramagnetic and AF phases is a fit to the data points given 
by Taf ^ 4:.5{x - Xc)^^'^, where Xc = 0.79. In the inset, 
versus x for T = 0.2. 
and 8 respectively. 



□ , o, and A, stand for L ■ 



nil 
24,20, 16, 



Let us explain now how to was extracted. To make 
sure that equilibrium was reached, plots of (72 and energy 
(average of H) were made over time intervals [t, 1.2t], 
not starting at t = to, as we do everywhere else, but 
starting at i = 0, from an initial random spin configu- 
ration. Semilog plots of 92 versus t displayed in Fig. [1] 
for X = 0.5, L = 24 and low temperatures show that 
a stationary state is reached only after some millions of 
MC sweeps. In order to check whether this state is truly 
in equilibrium, we define a time dependent spin overlap 
q, not among pairs of identical systems, but among spin 
configurations of the same system at two different times 
to and ti = to + t oi the same TMC run. 



q{to,t) ^ N-^Y.<Jj{to)(Jj{to+t). 



(6) 



Let q2{to,t) = {(q(to,t)) ). Suppose thermal equilib- 
rium is reached at a time t'. Then, 92(^01 ^ Q2 as 
t ~^ t' provided that to ^ t' ■ Plots of q2{to,t) vs t, for 
10~^to < t < to and to = 8 x 10^ MC sweeps, are shown 
in Fig. [1] Note that both quantities, 92 and $2 become 
approximately equal when t > 10^ MC sweeps. In or- 
der to obtain equilibrium results, we have always chosen 
sufficiently large values of to to make sure that, within 
errors, §2(^0: t) = Q2 for t^t^o- All values of to are given 
in Table I. 

In addition, we calculated a so-called correlation 
length for finite systems, 



1 



2sin(fc/2) 



92 



(I 9(k) 



- 1 



1/2 



(7) 



where 



k = llkil 



(8) 



Tj is the position of site j, and k = (27r/L, 0). Recall 
that this system is anisotropic, as interactions along the 
spin axis are twice as large as in the perpendicular direc- 
tion. Then, one could define a correlation length along 
the Y axis, £,y.Lj by choosing k = (0, 27r/L). We have 
found that is more convenient because it is less af- 
fected by finite-size effects than £,y,L- In order to com- 
pare with similar quantities defined for isotropic systems 
like the short-range 2D EA Ising model, we define also 
£,L — {£,x.L+£,y,L)/2. In contrast to P{q) and its first mo- 
ments, takes into account spatial variations of the EA 
overlap q and shows a good signature of SG transition. 
Its use has become customary in recent SG worki^i^ 
Analogous expressions define the AF correlation length 

f.{m) 



S,)^ by substituting for ijjj 



-1)^-(T, in Eqs.dUll). 



It is worth mentioning that in the ^l/L — > limit, 
is, up to a multiplicative constant, the spatial correlation 
length of (00 ^r)- Therefore in the paramagnetic phase 
we can think of , the L — > 00 limit of , as the true 
correlation length of a macroscopic system. On the con- 
trary, if there is strong long-range order with short-range 
order fluctuations (as predicted for the droplet model^i 
for 3D SGs), 92 7^ (that is, (0o0r) does not vanish as 
r — >■ 00) and {4>o4'r) — {4'o){4>r) would be short-range. It 
then follows from its definition Eq.(I7]) that ^ in 
2D. Following current usage, we shall nevertheless refer 
to as the "correlation length" . 
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FIG. 3: (Color online) (a) Squared staggered magnetization 
7712 vs T for X = 0.7. Icons o, □, o, and A stand for 
L = 24, 20, 16 and 8 respectively. Lines are only guides to 
the eye. Note that m2 decreases with L at all temperatures 
consistently with absence of AF order, (b) Same as in (a) 
but for X = 0.86. Note that m2 grows with L at low tem- 
perature, indicating an AF phase, (c) Same as in (a) but 
for the SG overlap parameter §2. (d) Same as in (c) but for 
X — 0.86. Direct comparison of curves shown in panels (b) 
and (d) for x — 0.86 indicate a coupling between m2 and q2. 
This coupling does not occur for x = 0.7 (see panels (a) and 
(c)). 




FIG. 4: (Color onhne) (a) Semilog plots of Cl/L versus T 
for x = 0.8, and L = 24 (o), L = 20 (□), L '= 16 (o), and 
L = 8 (A). In the inset, kurtosis of the m distribution versus 
T for the same values of x and system sizes, (b) Same as in 
(a) but for x = 0.86. 



III. RESULTS 
A. The AF phase 

The phase diagram shown in Fig. [5] summarizes our 
main results for the diluted 2D PAD model. For x > Xc 
we find a thermally driven second order transition be- 
tween the paramagnetic and AF phases at a Neel tem- 
perature T]^{x) that vanishes as x — >■ from above. 
The phase boundary meets the T = line at Xc — 0.79. 
For concentrations well below Xc the paramagnetic phase 
covers the whole range T >0. We do not find evidence of 
a SG phase at finite temperature. However, our results 
are consistent with a SG correlation length that diverges 
algebraically near or at Tsg = 0. In the following we 
report the numerical evidence on which these qualitative 
results are based. 

First we focus our attention on the paramagnetic- AF 
transition)^ The AF phase is defined by the staggered 
magnetization ([3]). We illustrate in Fig. [3|i how the mo- 
ment of staggered magnetization TO2 behaves with tem- 
perature for X = 0.7. Note that m2 appears to decrease 
as L increases even at low T. Plots of TO2 versus L (not 
shown) indicate a faster than algebraic decay in L, as one 
expects for a non AF phase. This is in sharp contrast to 
the behavior of ?7i2 for x ~ 0.86 (see Fig. [3b). Curves for 
different L cross at Tn ~ 1.15. Below this temperature 
m2 increases with L indicating the existence of an or- 
dered AF phase. Similar results are obtained for higher 
values of x. In the inset of Fig. [21 plots of mi versus x 
for different system sizes at low temperature show that 
the system exhibits AF order for x > 0.8. Similar graphs 
were obtained for m2- These are our first pieces of evi- 
dence for the existence of an AF phase above Xc ^ 0.8. 
It is instructive to compare the behavior of m2 with that 
of q2 shown for x = 0.86 in Figs. [8)3 and [Sji. TO2 and (72 
are not qualitatively different. This is not so for a; = 0.7 
where there is no AF order (compare Figs. [3^ and [3}:). 
From Fig. [3j; it is not obvious whether 52 vanishes or not 
as L increases at very low temperatures. We will return 
to this point in the discussion of Fig. [5l 

For further information about the extent of the AF 
phase, we also examine how the cumulant-like quantity 
gm and the finite-size AF correlation length behave for 
several pairs of values x, T. Let us first outline how 
is expected to behave in the various magnetic phases. It 
clearly follows from its definition that g^. — > 1 as i 
00 in the case of long-range AF order. From the law 
of large numbers it also follows that -> as L — ?> 
00 in the paramagnetic phase. These two statements 
imply that curves of g^n vs T for various values of L cross 
at the phase boundary between the paramagnetic and 
AF phases. We make use of this fact to quantitatively 
determine the paramagnetic-AF phase boundary. Plots 
of g,n vs T are shown in the insets of Figs. [4^ and[4|D for 
X = 0.8 and 0.86, respectively. The signature of an AF 
phase below T ~ 1.2 is clear for x = 0.86. The inset of 
Fig. [3^ shows that within errors curves of gm vs T for x = 
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0.8 and different system sizes merge instead of crossing 
at and below T = 0.5(1). Note also that Qm does not go 
to 1 as T — 7> 0, indicating a broad distribution of P{m) 
even in this limit. Finally, for x < Xc (see Section UlI Bp . 
we find that decreases as L increases for all T which 
is consistent with the absence of AF in this region. 

In recent literature on SG phases the scale invariant 
finite-size correlation length^ is frequently used to give 
evidence for a finite temperature transition, since S,l/L 
crosses at the transition temperature Tjv and spreads out 
above and below T/v. The advantange of ^l/L over kur- 
tosis is that the former may even diverge as L — > oo in 
contrast to the latter that tends to 1. Then we use the 
AF correlation length ^^'"^/L to pinpoint values for Ti 



N 



by the value of T where curves cross. Recall that <^j^'"'' be- 
comes a true correlation length when ^j™^ /L ^ 1 . Then, 
in the paramagnetic phase, 'Cl"V^ 0(11 L)^ therefore 
decreasing as L increases. At T/v, S^l^^ /L must become 
size independent, as expected for a scale-free quantity. 
At lower temperatures, well in the long-range AF phase, 

we expect ^["V-^ ~ 0{L). Plots of Ci"V-^ versus T 
are shown in Figs. and [4|3 for x = 0.8 and 0.86, re- 
spectively. Note that curves spread out above and below 
Tn = 1.20(5) for x = 0.86. Similar graphs for x = 0.9 
allow one to obtain the value T/v = 1.50(5). 

On the other hand curves merge for all temperatures 
below T = 0.5(1) for x ^ 0.8 and i > 16 (sec Fig. Ha,), 
while TO2 decreases, within errors, algebraically with L for 




FIG. 5: (Color online) (a) Plots of q2 versus L for x — 
0.2. o, □, o, l>. A, V, o, •, and x stand for T/x — 
0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, and 1.0 respectively. Lines 
are to guide the eye. Clearly, data for T/x > 0.4 deviate from 
the straight dashed lines implying faster than a power of 1 /L 
decay, (b) Same as in (a) but for x = 0.5. Here, only data for 
T/x > 0.7 decay faster than a power of 1/L. For all data, we 
have checked that, within errors, Q2 = q2- 



the studied system sizes (not shown). It is interesting to 
note that graphs of the SG quantities g and ^l/L (instead 
of gjn and the AF S^l^^ /L) give qualitatively the same 
picture when plotted versus T except from the fact that 
g — >■ 1 as T — >• for all x. Thus, the most straightforward 
interpretation of the data shown in Fig. [4^ is that for all 
temperatures below T = 0.5 the system is near or at the 
AF phase boundary and its behavior displays criticality. 

We have thus established all points of the AF phase 
boundary shown in Fig. [5]for x > 0.7. A fit to these data 
points, given by T/v — 4.5(a; — 2;^)^^^, where x^ = 0.79(5) 
is shown in Fig. O Finally, for x < 0.7 (see below) we 

find that ^j^"'''/T decreases as L increases for all T, as 
expected. 



B. Very diluted systems 

This Section is devoted to show the numerical results 
drawn for distributions of q and their first moments, and 
for for systems with weak concentration. As for 3D 
PADs,— we expect universal behavior for a; ^ 1 which 
enables us to compare our results with previous work. 
Thus, we direct our attention on the results we have ob- 
tained for x = 0.2 and x = 0.5. Both values are well 




-0.5 0.5 

q, m 



1 1 
q',m' 



FIG. 6: (Color online) (a) Plots of the probability distribu- 
tions P{q) versus q, and P{m) versus m both for x — 0.2 and 
T/x = 0.4. o, □, o are for P{q) and system sizes L = 32,24 
and 20, respectively. A, V, and x are for P{m) and system 
sizes L — 24, 20 and 16, respectively, (b) Same as in (a) but 
for the scaled distributions P{q') versus q' = q/qi, and P{m') 
versus m' = m/m\. The thick dashed line corresponds to 
the Gaussian distribution of paramagnets in the macroscopic 
limit, (c) Same as in (a) but for x = 0.5 and T/x = 0.4. o, 
□ , o are for P{q) and system sizes L — 24, 20 and 16, respec- 
tively. A, V, and X are for P{q) and system sizes L = 24,20 
and 16, respectively, (d) Same as in (b) but for x = 0.5 and 
T/x = 0.4. 
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FIG. 7: (Color online) Kurtosis g of the distribution of q ver- 
sus T/x for X = 0.2. •, o, □, A, and o stand for system sizes 
L = 32, 24, 20, 16, and 8 respectively. In the inset, semilog 
plot of the relative mean square deviation versus T/x for 
the same x and system sizes. Lines are guides to the eye. (b) 
Same as in (a) but for x = 0.5. o, □, A, and o, stand for 
L = 24, 20, 16, and 8 respectively. In the inset, scaling plot 
{g as a function of {T j x)L}^ of the data shown in the main 
figure. 

below the x> Xc — 0.79 region, in which AF appears at 
low temperatures. 

A plot of versus T for x = 0.7 is shown in Fig. [Sj;. 
Qualitatively similar graphs are obtained for other values 
of X satisfying x Xc- Note that qi decreases as L in- 
creases, even at low temperatures. It is difficult to deduce 
from these plots whether or not qi vanishes as i — > oo 
at low T. In order to elucidate this question we prepare 
log-log plots of 52 vs L, shown in Figs. [5^ and for 
X = 0.2 and 0.5 respectively. Data points in these figures 
seem to be consistent with a decay faster than (72 ^ 
for T/x > 0.3, indicating that we are in the paramag- 
netic phase. Plots of qi vs L show the same qualitative 
behavior. Altogether these results leave small room for 
the existence of a SG phase with quasi-long range order 
at very low temperatures, as it has been reported for the 
3D diluted PAD model for T/x > 1. 

Next we report results for the distributions of m and 
q at low temperature. Due to the central limit theorem 
and since the correlation lengths are finite in the param- 
agnetic phase, P{m) and P{q) are expected to be normal 
distributions for L = 00. The droplet picture^! for SG's 
predicts that P{q) = [5{q -\- qo) + S{q — 9o)]/2 where go is 
the EA order parameter, and that the tail of P{q) down 
to g = for finite-size systems vanishes as L increases. 
On the contrary, the RSB picture^ predicts a nontriv- 




0.5 1 1.5 2 

T/x 

FIG. 8: (Color online) (a) Semilog plots of (a) SG correlation 
length divided by system size £,x.l/L versus T/x for x — 0.2, 
and L = 32 (•), L = 24 (o), L = 20 (□), L = 16 (A), L = 12 
(o) , and L = 8 ( X ) . Lines are guides to the eye. (b) Same as in 
(a) but for x = 0.5, and L = 24 (o), L = 20 (□), L = 16 (A), 
and L = 8 (o). In the inset of Fig. (a), log-log plot of ^x,l/L 
versus 1/L for x = 0.2. o, □, o. A, •, <, x, and ■ stand for 
T/x = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 1.0 respectively. 



ial distribution with a non- vanishing P{q = 0) which is 
size independent. Plots of P{q) and P{m) are shown for 
X = 0.2 and x = 0.5 for T/x = 0.4 in Figs. EK-Et. All 
distributions depend on L. Pirn) are found to be nor- 
mally distributed. In Figs. [6)3 and[6jl normalized distri- 
butions of the reduced quantity m' = m/mi are shown. 
Note that P{m') ~ (I/tt) exp(— to'^/tt) for the studied 
system sizes, indicating complete absence of AF order. 
On the other hand, P{q) are found to be double peaked 
distributions. As L increases, their peak positions shift 
towards q = Q and P(0) increases. Neither the droplet 
nor the RSB models consent a fit to these data. If it 
turns out that our systems are near or at criticality, then 
P{q') (where g' = q/qi) ought to be size independent. 
However, reduced distributions P{q') m Figs. [n}3-[5Ji are 
shown to have an L dependence. We conclude that our 
results are only consistent with a paramagnetic phase. 
Similar conclusions apply for T/x > 0.2. 

In the same way as explained in Section Fill Al for quan- 
tities gm and ^^'"■'/L, their dimensionless SG counter- 
parts g and S,l/L, as well as A^, indicate the location 
of the temperature Tsg of a SG transition. Recall that, 
according to the finite-size scaling assumption, all these 
quantities depend only on L{T ^TsgY then become 
size independent at Tsg- We are assuming that for large 
enough sizes, L/S^ (where ^ is the true correlation length) 
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is the only relevant parameter and £, ^ {T ~ Tsg)~'^ ■ 
Plots of g vsT/x are given in Fig. [TK and[7|3 for x = 0.2 
and 0.5 respectively. It seems that curves of g for various 
values of L do not cross and only merge as T — > 0. This is 
consistent with Tsg = in accordance with the behavior 
of 2D EA systems, although a merging at T/x « 0.2 is 
not completely excluded within errors. We found more 
useful to study A^, which has a direct interpretation 
as the uncertainty of q/qi and could be computed with 
higher precision as it involves lower moments P{q)- 
A plot of vs T/x for x = 0.2 is shown in the inset of 
Fig. [7^. Clearly, A^ increases with L for all T/x> 0.2 as 
expected for a paramagnetic phase and contrary to the 
expected behavior A^ — >■ for increasing L in a SG phase 
with non- vanishing order parameter. A similar behavior 
for X = 0.5 is observed (not shown). 

We now turn our interest to the behavior of the SG 
finite-size correlation length £,l/L. Similarly to g, S,l/L 
is independent of L at Tsg as it corresponds to be for a 
scale-free quantity. On the other hand, well inside the 
paramagnetic phase, ^l/L should diminish as 0{1/L), 
provided £,l/L ^ 1. Figs. [8^ and [SJd show data for 
S,x.L as a function of T/x for several system sizes, for 
X = 0.2 and 0.5 respectively. Curves do not cross at any 
finite temperature and £.x,l/L decrease as L increases at 
least for T/x > 0.4. All curves seem to converge only 
as T — > suggesting Tsg = 0. Plots for ^x,l vs L (see 
inset of Fig. [8^) are consistent with an algebraic decay 
£,x,l/L ^ L~y for temperatures T/x > 0.4 and system 
sizes L > 16. For lower temperatures, data do not vary 
very much as L increases and so we cannot definitely 
rule out a non-vanishing ^x.l/L in the thermodynamic 
limit. Consequently, from Fig. [5] alone, a transition at 
very low, but not zero, temperature cannot be completely 
excluded. 



C. The u exponent at Tsg = 

Our results concur with the behavior found 
for a random-bond 2D RAD model with dipolar 
interactions,— and the one found time ago for 2D EA 
modelsj^ Recent simulations for the latter (with larger 
system sizes and lower T) including nearest-neighbor 
exchange interactions find a lower value 1/v = 0.29(4) 
for Gaussian distributions^^ but do not provide con- 
clusive results for bimodally distributed interactions 
A scenario has been proposed where these models for 
varying realizations of disorder belong to the same 
universality class at non-zero temperatures?^ even 
though their respective behaviors differ at T = 0. 

Let us assume that Tsg = and that ^ diverges as 
^ ^ T~'^ . According to finite-size scaling, dimcnsionless 
quantities like g and / L should scale as 
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FIG. 9: (Color online) (a) Kurtosis y as a function of the 
finite-size correlation length divided by system size, ^l/L for 
X = 0.2. •, o, □, o, and A stand for L = 32, 24, 20, 16, and 12 
respectively. All data should collapse onto an universal curve, 
provided that scaling corrections are small. The thick contin- 
uous line stands for the universal curve that corresponds to 
the 2D Ising SG model with short-range interactions, (b) 
Same as in (a) but for x = 0.5. o, □, o, and A stand for 
system sizes L = 24, 20, 16, and 8 respectively. 
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g = G{TL^/'^), iL/L = X{TL'/'^) 
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FIG. 10: (Color online) (a) Scaling plot for the SG correlation 
length divided by system size, (,x,l/L versus {T/x)L'^''' with 
l/v = 0.35 for X = 0.2. •, o, □, and A are for L = 32, 24, 20 
and 16, respectively, (b) Same as in (a) but for x = 0.5. o, 
□ , A, and o are for L = 24, 20, 16 and 8, respectively. Error 
bars, where not shown, are smaller than symbols. 



It follows that g = F{^i^/L), where F is in principle 
a non-universal function that, apart from the bulk uni- 
versality class, depends on the boundary conditions (we 
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chose them to be periodic), the sample shapes (squared 
lattices), as well as the anisotropy^ of the interactions. 
Plots of g versus £,l/L are shown in Figs. and [5)3 for 
X = 0.2 and 0.5 respectively. Data for different values 
of L and x should collapse into a single scaling curve, 
on condition that finite-size effects are small. This is 
what we find for systems with > 100. However, for 
smaller N curves spread out indicating that finite-size 
scaling corrections are large. It is remarkable that data 
seems to collapse onto the universal scaling curve found 
for the (isotropic) 2D EA models mentioned above in this 
Section^i2^ specially at low T. Note that in Figs. [SJi 
and|9)D we use £,l/L instead of ^x,l/L in order to aver- 
age anisotropic effects over the two principal axes of the 
underlying square lattice. This suggests that both, 2D 
PAD and short-range EA models, may share a common 
universality class. However, large correlations to scaling 
for the available system sizes prevent us to go further on 
this direction. 

Scaling plots of S,x.l/L versus [T /x)L^/^ are shown in 
Figs. [TOk andfTOb for x = 0.2 and x = 0.5 respectively. 
Due to the presence of large finite-size corrections, no 
value of 1 allows to collapse all data in one single curve. 
We have chosen l/i/ in order to allow data collapse for 
large L and low T. We find that = 0.35 is a suitable 
value in both cases. Scaling plots become significantly 
worse (not shown) when using values of l/v outside the 
interval [0.25,0.45] ([0.3,0.4]) for x = 0.2 {x = 0.5). That 
gives a rough estimate of the error on 1 /u. A scaling plot 
of g versus {T/x)L^^'^ is shown in the inset of Fig. [Tb for 
X ~ 0.5. Again, l/iy = 0.35(10) allows to scale data for 
large L and low T, which is consistent with the value 
extracted from ^x,l/L. The value l/v = 0.35(10) also 
agrees with the effective exponent found in old simula- 
tions of the 2D EA model for small system sizes and 
relatively high temperatures and is slightly larger than 
the value l/iy = 0.29(4) found in recent simulations for 
the same models (but the two values are still consistent 
within errors). In any case, such discrepancies may be 
caused by the fact that finite-size scaling corrections arc 



important for the limited system sizes that we were able 
to simulate. 



IV. CONCLUSIONS 

By tempered Monte Carlo calculations, we have stud- 
ied a diluted dipolar Ising model on a square lattice. 
There are only dipolc-dipolc interactions. Spins ran- 
domly occupy only a fraction x of all lattice sites. The en- 
tire phase diagram of the system, in particular the bound- 
ary between the AF and the paramagnetic phases, has 
been explored and it is shown in Fig.O We have also pro- 
vided strong evidence that the paramagnetic phase covers 
the whole T > range for x < Xc, where Xc = 0.79(5). 
From the behavior of the spin-glass (SG) overlap q, the 
relative mean square deviation A^, kurtosis g, and £,l/L, 
we conclude that Tsg = for a; < a;c, and there is an al- 
gebraic divergence of the correlation length with an expo- 
nent l/i' = 0.35(10). All these properties are consistent 
with the behavior found for the 2D diluted RAD and 2D 
EA model with short-range interactions. This is to be 
contrasted with the manifestly different behavior found 
for the 3D PAD (quasi long-range order at low T) and 
the 3D EA model (with a non-vanishing order parameter 
in the SG phase). 
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